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ABSTRACT 

We treat the problem of adiabatic losses and stochastic particle acceleration in gamma-ray 
burst (GRB) blast waves that decelerate by sweeping up matter from an external medium. The 
shocked fluid is assumed to be represented by a homogeneous expanding shell. The energy lost 
by nonthermal particles through adiabatic expansion is converted to the bulk kinetic energy of 
the outflow, permitting the evolution of the bulk Lorentz factor T of the blast wave to be self- 
consistently calculated. The behavior of the system is shown to reproduce the hydrodynamic 
self-similar solutions in the relativistic and nonrelativistic limits, and the formalism is applicable 
to scenarios that are intermediate between the adiabatic and fully radiative regimes. 

Nonthermal particle energization through stochastic gyroresonant acceleration with magnetic 
turbulence in the blast wave is treated by employing energy-gain rates and diffusive escape 
timescales based upon expressions derived in the quasilinear regime. If the magnetic field in the 
shocked fluid approaches its equipartition value, this process can accelerate escaping particles to 
> 10 20 eV energies, consistent with the hypothesis that ultra-high energy cosmic rays (UHECRs) 
are accelerated by GRB blast waves. Due to particle trapping by the magnetic turbulence, only 
the highest energy particles can escape during the prompt and afterglow phases of a GRB for 
acceleration by a Kolmogorov spectrum of MHD turbulence. Lower energy particles begin to 
escape as the blast wave becomes nonrelativistic and shock Fermi acceleration becomes more 
important. 

Subject headings: acceleration of particles — cosmic rays — gamma rays: burst 



1. Introduction 

The origin of cosmic rays (CRs) with energies above the knee of the CR spectrum at 3 x 10 15 
eV is not established. Nor is the origin of ultra-high (> 10 19 eV) energy cosmic rays (UHECRs) known. 
Because of their large Larmor radii, UHECRs probably originate from extragalactic sources. Theoretical 
difficulties (Lagage and Cesarsky 1983) to accelerate cosmic rays to energies above the knee energy via shock 
Fermi acceleration in supernova remnants suggests that a new class of sources with adequate power and 
acceleration efficiency is required. The stellar progenitors of gamma-ray bursts (GRBs) provide a plausible 
solution (Milgrom and Usov 1995, 1996; Waxman 1995; Vietri 1995). One motivation for this proposal is the 
coincidence (Waxman 1995; Vietri 1995) between the power supplied by GRB sources within the Greisen- 
Zatsepin-Kuzmin photopion energy-loss radius and the power needed to produce the observed energy density 
of UHECRs. This coincidence has recently been verified (Dermer 2000) for the external shock model of GRBs 
(Bottcher and Dermer 2000), under the assumption that the rate density of GRBs follows the star formation 
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history of the universe. A second motivation is that relativistic shock waves provide a site to accelerate 
particles to much higher energies than is possible in nonrelativistic shock waves. For example, reflection 
from a relativistic shock with Lorentz factor T will increase a particle's energy by a factor ~ T 2 . GRB blast 
waves with 100 < T < 1000 can therefore accelerate particles to energies near the knee of the cosmic ray 
spectrum in a single cycle, with subsequent acceleration producing cosmic rays above the knee energy (Vietri 
1995; Milgrom and Usov 1996). 

GRB emission originates either from internal shocks in a relativistic wind (Rees and Meszaros 1994) or 
from an external shock (Rees and Meszaros 1992; Meszaros and Rees 1993) that forms when a relativistic 
blast wave decelerates and is energized by sweeping up matter from the surrounding medium. Acceleration 
to energies above the knee of the CR spectrum depends on the specific acceleration mechanism. Gallant and 
Achterberg (1999) show that an external shock model to accelerate UHECRs starting from nonrelativistic 
particles is not viable. In this paper, we instead consider a scenario for UHECR acceleration based upon 
stochastic particle acceleration in GRB blast waves that slow from relativistic to nonrelativistic speeds (Wax- 
man 1995; Rachen and Mesazaros 1998; Schlickeiser and Dermcr 2000). Although we focus on gyroresonant 
acceleration in relativistic shocks formed by the sources of GRBs, this approach is also applicable to stochas- 
tic acceleration in systems with mildly relativistic and nonrelativistic shocks. For example, Type Ib/c SNe 
produce shocks with speeds that often exceed 0.1c (Weiler et al. 2000). Mildly relativistic shocks with Lorentz 
factor r ~ 1.6-2 were produced by the Type Ic SN 1998bw (Kulkarni et al. 1998). Irrespective of whether 
GRB 980425 is associated with SN 1998bw (Pian et al. 1999), the existence of SN 1998bw demonstrates 
that unusual SNe such as SN 1998bw eject relativistic flows that will contribute to CR particle acceleration. 
The adiabatic and stochastic processes considered here will also play a role in particle acceleration in SNRs, 
although the first-order systematic energy-gain rate dominates that of stochastic gyroresonant acceleration 
in nonrelativistic shocks. 

In Section 2, we treat adiabatic losses in GRB blast waves, and provide an equation for blast-wave 
evolution that accounts for these losses in a manner that self-consistently rechannels the energy lost from 
particle expansion into the directed kinetic energy of the expanding blast wave. An accurate treatment 
of adiabatic losses is required to estimate maximum particle energies through stochastic acceleration. We 
assume throughout that the blast wave is uniform and can be described by a thin shell of outflowing matter, 
whereas in reality the flow is subject to instabilities (Kang, Jones, and Ryu 1992) and nonlinear feedback from 
the accelerated particle distribution (e.g.,Blandford and Eichler 1987; Baring et al. 1999) that can smooth 
the transition at the shock front. Acoustic instabilities can also cause a turbulent fragmentation of the 
shocked fluid (Ryu and Vishniac 1987). Although the turbulence so generated can be effective in generating 
a magnetic field near equipartition with the downstream flow, it also invalidates the simple assumption of 
a uniform shocked shell, thus limiting the applicability of the method. Future work must consider how 
instabilities and nonlinear effects of the shock structure will modify the results presented here. 

Stochastic particle acceleration through gyroresonant acceleration is considered in Section 3 by general- 
izing quasilinear expressions to a fully turbulent regime with relativistic Alfven speeds, as probably applies 
to GRB blast waves. We show that UHECR acceleration is possible in relativistic blast waves by considering 
various limits to particle acceleration, including a comparison of the particle gyroradius with the blast wave 
width (Hillas 1984), and competition with adiabatic losses, synchrotron losses, and diffusive escape. Our 
work improves upon previous treatments that have made similar comparisons (Vietri 1995, 1998a; Rachen 
and Mesazaros 1998), though here we restrict our considerations to stochastic particle acceleration in GRB 
blast waves. Particle energy evolution and escape from the blast wave is examined through the use of a 
particle continuity equation. 
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Simulations of the spectra of particles accelerated in GRB blast waves are presented in Section 4. 
Particles accelerated through stochastic processes in relativistic blast waves will also participate in shock 
Fermi processes as the blast wave decelerates to nonrelativistic speeds, so that a full calculation of the 
spectrum of the accelerated particles must consider both first- and second-order Fermi acceleration. GRB 
sources may thus be the sources of UHECRs and CRs above the knee of the CR spectrum, as well as 
contributing some portion of lower energy CR hadrons. Summary and conclusions are given in Section 5. 



2. Evolution of Blast Wave Lorentz Factor 

Both the particles that are captured from an external medium by a GRB blast wave and the ejected 
thermal plasma particles that initially carry most of the explosion energy will be subject to adiabatic losses 
due to volume expansion of the blast wave shell. This energy is reconverted into the directed kinetic energy 
of the outflow. Here we present a treatment of expansion losses in a blast wave that contains particles with 
arbitrary energies in the comoving frame. The treatment is valid for a blast wave with general values of the 
bulk speed, and is therefore applicable to both decelerating relativistic GRB blast waves and to nonrelativistic 
SNR shock waves. The treatment applies to general deceleration regimes from the fully adiabatic to the fully 
radiative limit. 



2.1. Blast- Wave Equation of Motion 

Let p' = \/"f 12 — 1 represent the dimensionless momentum of a particle with Lorentz factor 7', where 
primes refer to quantities in the comoving frame. The particle distribution function N'(p';x), integrated 
over the volume of the blast wave shell, is defined so that N'(p';x)dp' — [dN'(p';x)/dp']dp' represents the 
differential number of particles with momenta between p' and p' + dp' in a blast wave at radius x. In writing 
this function, we assume that the particles are isotropically distributed in the comoving frame. Global 
conservation of energy implies that 

d{TU t ot) = dm + TdU rad . (1) 

All masses are in energy units. Here T is the Lorentz factor of the blast wave, dm is the differential change 
in the swept-up rest-mass energy, Utot is the total internal energy in the comoving frame, and dU ra d is 
the differential change in the comoving internal energy that is radiated from the blast wave. The radiated 
energy is assumed to be isotropically emitted in the comoving frame. We do not take into account energy 
gains when the system scatters or captures external radiation or magnetic-field energy. Thus equation (1) 
can apply to internal emission processes such as synchrotron, synchrotron self-Compton, brcmsstrahlung, 
photomeson and secondary nuclear production processes, but not to Compton scattering processes involving 
external photon fields. 

The quantity 

Utot = Hm p / dp'iN\p'- x) = M + U = M + m(x) + U (2) 
Jo 

represents the total internal energy, including rest-mass energy, whereas 

/>oo 

U = i im v \ dp\i -l)N'{p';x) (3) 
Jo 

represents the internal energy with rest-mass energy excluded. The total rest mass M consists of the sum of 
the mass M = E /T ejected by the explosive event with energy E and initial Lorentz factor r , in addition 
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to the swept-up mass m(x) — f£ dx [dm(x)/dx] — A-K^m p J Q X dxx 2 n ext {x). Here n ext {x) is the density of the 
external medium, and is assumed for simplicity to have radial symmetry. The quantity fi is the ratio of the 
mass of swept-up particles to the mass of swept-up protons, assuming that the different ionic species have 
the same comoving distribution function. In most case, the protons make the dominant contribution to the 
swept-up nonthermal particle energy; thus /i = 1. Generalization to multiple species of swept-up particles, 
including leptons, ions and charged dust (Schlickeiser and Dermer 2000) is straightforward but is not treated 
here. The explosion is assumed to be uncollimated, though it is also straightforward to generalize to a system 
with beamed energy releases and to external media with arbitrary density distributions. 

The internal energy U in the blast wave increases by the addition of the kinetic energy of the swept- 
up matter, and decreases due to particle energy losses through adiabatic expansion and radiation. Thus 
dll = dllrn + dUadi + dU ra d, where dU m = (r — l)dm and dU a di represent the changes in internal energy 
due to swept-up particle kinetic energy and adiabatic losses, respectively. Expanding equation (1) gives 
d[T(U + M)] = TdU + Tdm + (U + M)dT. From this follows the equation of blast wave evolution 

dT _ dm + (^2)dU a di 



r 2 - 1 m + u 



where the blast wave momentum P = \/T 2 — 1. Equation (4) can equivalently be written as 

_ dP = PT(dm/dx) + (^){dU adt /dx) 
dx ~ M + m{x) + U 



(5) 



The second term in the numerator on the right-hand side of either equations (4) or (5) represents the 
impulse to the blast wave resulting from adiabatic energy losses of the particles in the blast wave. Equation 
(4) generalizes the equation derived by Panaitescu et al. (1998) for an adiabatic blast wave by the inclusion 
of radiative losses in the calculation of U and dU a di/dx. 

For a strongly radiative blast wave, the adiabatic loss term is small and the radiative equation of blast- 
wave evolution derived by Blandford and McKee (1976) is recovered by setting U = dU a( n/dx = and 
dm = dM, giving 

dT dm , . 

(6) 



T 2 - 1 M + m(x) 



This is easily solved to obtain 



r( , = [i + m(z)/M ] 2 (r + i) + r -i 

1 ' [l + m(a;)/Mo] 2 (ro + l)-ro + l ' 1 ' 

When radiation losses are included but adiabatic losses are neglected, equation (4) reduces to the form 

dT dm 



T 2 - 1 M + m(x) + U 



(8) 



that is derived from a momentum- and energy-conservation analysis (Dermer and Chiang 1998; Chiang 
and Dermer 1999; Piran 1999). This equation can be solved (Chiang and Dermer 1999) by noting that 
dll = (V— \)dm. Furthermore, if a fraction e of the swept-up kinetic energy is instantaneously radiated from 
the blast wave, then dU = (l-e)(T-l)dm. Defining Al = M a +m(x) + U implies that dM = [e+T(l-e)]dm, 
one obtains 

dT _ dM 
~ 1 ~ M[e + T(l-e)] ' { ' 
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Equation (9) has been used to treat partially radiative blast waves (Piran 1999; Bottcher and Dernier 
2000; Moderski et al. 2000). Huang, Dai, and Lu (1999, 2000) suggest using the relation dM = dm[2r(l-e) + 
e] so that the blast- wave evolutionary behavior also follows the Sedov behavior in the adiabatic nonrelativistic 
regime. Both approaches are adequate to examine the approximate behavior of relativistic blast waves in 
intermediate radiative regimes, and the treatment of Huang and colleagues furthermore approximates the 
correct dependence of nonrelativistic adiabatic blast waves. But these approaches do not properly treat 
adiabatic losses of particles within the blast wave. Moreover, it is necessary to make the assumption that 
e w const throughout the regime of blast-wave evolution under consideration. When questions of particle 
acceleration are treated, it is essential to have a correct treatment of adiabatic losses, which equation (5) 
provides. 



2.2. Adiabatic Particle Losses 

If the internal energy U in the blast wave is dominated by the kinetic energy of a thermal particle 
distribution, then U changes due to volume expansion according to the relation U~ l {dU a di/ dx) — — (7 — 
l){d\nV' / dx), where 7 is the ratio of specific heats, and 7 = 4/3 and 5/3 for relativistic and nonrelativistic 
monatomic gases, respectively, and V is the comoving fluid volume. The general case involving a mixed 
fluid that includes both thermal and nonthermal particles can be derived by noting that p' changes through 
expansion losses according to the expression 

dp' _p'd\nV 

- { ^ )adi -^^r- (10) 

Equation (10) reproduces the limiting adiabatic loss forms exhibited by monatomic thermal gases, noting 
that in the relativistic limit p' » 1, dry 1 — > — (7'/3)dln V 7 , and in the nonrelativistic limit p'<l, d(p' 2 /2) — > 
-(2/3)(p' 2 /2)dlnV. 

The comoving volume of the fluid shell is approximated by the expression 

V = Anx 2 A' =4ttx 2 (^) , (11) 

where A' is the comoving width of the blast-wave shell. As noted in the Introduction, this is an extreme 
simplification due to instabilities and turbulence that will modify the shell structure (Ryu and Vishniac 
1987; Meszaros, Rees, and Papathanassiou 1994). When /a = 1/12, this relation is consistent with the 
approximation V = m(x)/(n'm p ) (Panaitescu et al. 1998) in the limits T » 1 and T — 1 <C 1, where 
n'(x) = (7T + l)n ext /(j — 1) is the downstream comoving density in terms of the external medium density 
n ext . The differential adiabatic expansion energy- loss rate for particles in a GRB blast wave is therefore 
given by 

.dp'. ,,\ ldlnT. 

-f )adi = -p' ----7- > 12 

dx x 6 dx 

using equations (10) and (11). Note that the second term on the right-hand-side in the parentheses of 
equation (12) is small in comparison with the first term in the limit r — 1 <C 1. 

The internal energy therefore changes through adiabatic losses according to the relation 
Modi A ldlnT [°° , p' 2 

— = - m *> i x-3— ) J, ^ ( V )N(p;x), (13) 
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letting = 1. If only adiabatic losses are important, then equation (12) is easily solved to give 

pi=p\x,x t )=pi(^){^}V\ (14) 

where p\ is the momentum of a particle injected at radius Xi when the blast wave was moving with Lorentz 
factor T(xi). 

The sweep-up function 

QM,Xi) = dN ^' l d ^ = 4 7 rx J 2 rw(x J ) ( 5[f4 - P{ Xi )\ (15) 

is defined so that Q su (p' i , Xi)dp' { dxi gives the differential number of particles with comoving- frame momentum 
p\ between p\ and p' i ^ r dp' i that are swept-up between radii a;, and Xi + dxi. The comoving particle distribution 
function at location x is therefore 



V (>':.»■)= / d.r H \^\Q su {p', Xi ). (16) 



f 

Jo 

From the definition of U in equation (3), we obtain 

U(x) — 47rm p n / dxixj (7 — 1) , (17) 
Jo 



where p = (xi/x)[T(x)/T(xi)] 1 ^ 3 P(xi) and 7 = \Jp 2 + 1. Here the 5-function in equation (15) is used to 
solve the integral over dp' in equation (3), and we simplify to a uniform surrounding medium with density 
no. Similarly, the dU a di/dx term in equation (13) is evaluated to obtain 

dU adl . ,\ ldlnr 



dx 



a 1 dinr\ r , 2 ,p 2 x 

47rm p n ( -— — ) / dx i x i ( — ) . (18) 

X o G(X ,/o 7 



Fig. 1 shows a calculation of equation (5) under different assumptions for the particle energy losses. 
The thermal ejecta particles are assumed to be cold, so that only energy losses of the swept-up particles are 
considered. The parameters of the calculation are E n = 10 54 i?54 ergs, r = 300r 30 o, and n — 100n 2 cm~ 3 , 
with £54 = r 30 o = ri2 = 1. For these parameters, the deceleration radius (Meszaros and Rees 1993) 

X d EE ( 3 f° )V3 S 2. 6 x 10 16 ( ^_)l/3 cm . (19) 

We note that models of afterglow spectra (Wijers and Galama 1999; Panaitescu and Kumar 2001) indicate 
that GRBs may take place in relatively tenuous environments with 10~ 4 cm~ 3 < no <^ 10 2 cm -3 , so 
that our standard parameters are on the high range of inferred values of no- These models have neglected, 
however, to consider adiabatic losses of the injected electron distributions, which are important for electrons 
with Lorentz factors that produce synchrotron emission near the cooling break. Moreover, the blast-wave 
dynamics do not take into account evolution in intermediate radiative regimes that would be important 
if UHECRs carry internal energy away from the system. Comparisons (Dermer et al. 2000) of detailed 
numerical models with analytic equations for afterglow emission show discrepancies by well over an order 
of magnitude, particularly when Compton losses are important, as is the case in many of the afterglow fits. 
It is not clear how these effects may change inferences about no, but the reader should keep in mind the 
possibility that n ~ 0.01-1 cm~ 3 is more realistic, and we will consider such densities in the subsequent 
discussion. 
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The numerical solution of equation (5), using equations (17) and (18) to calculate the various terms 
in this equation, is shown by the thick solid curve in the figure. It exhibits the behavior P(x) oc x~ 3 / 2 at 
x » Xd- This gives the T(x) oc x~ 3 / 2 dependence of T at relativistic speeds, and the Sedov solution behavior 
f3(x) oc x~ 3 / 2 at nonrelativistic speeds, where P — (3T. For comparison, the dashed lines show the asymptotes 
for the shocked fluids at relativistic and nonrelativistic bulk Lorentz factors (Blandford and McKee 1976). 
In the relativistic case, the Lorentz factor of the shocked fluid is given by T = (17E/16tt Po ) 1/2 x- 3 / 2 when 
x ^> Xd and r>l, where p is the rest-mass energy density of the external medium. In the nonrelativistic 
case, (3 = 0.29(E/ p^Y^x^ 3 / 2 , where we assume that 7 = 5/3 is the adiabatic index of the external unshocked 
gas, and we note that the speed of the shocked fluid is 3/4 of the speed of a nonrelativistic shock. 

The analytic solution to the equation of blast-wave evolution in the radiative regime (Blandford and 
McKee 1976), equation (7), is shown by the labeled solid curve for this case. When x » Xd, the P oc x~ 3 
behavior is recovered. The analytic solution to the equation for blast-wave evolution with momentum- 
conservation, equation (9), was derived by Chiang and Dermer (1999) and is shown by the second labeled 
solid curve in Fig. 1. When x 3> Xd, P oc x~ 3 l 2 when F>1 and P oc x~ 3 when P <C 1. As pointed out by 
Huang, Dai, and Lu (1999), this equation does not reproduce the Sedov solution at nonrelativistic Lorentz 
factors, but rather follows the momentum-conservation equation derived by Oort (Lozinskaya 1992). 



2.3. Evolution of Comoving Particle Distribution Function and Internal Energy 

When only adiabatic losses are important, the comoving particle distribution function is given by 

N'(p'; x) = 4irn [ d Xi x 2 D S[Dp' - P(xi)} , (20) 
Jo 

where D = (x/xi)[T(xi) /T(x)] 1 / 3 , using equations (14)-(16). The internal energy U a di, for a blast wave 
subject to adiabatic losses only, is given by substituting equation (20) into equation (3). 

The numerical results show that an adiabatic blast wave evolves according to the relation P(x) w P 
for x <C Xd, and P(x) w Po(x/xd)~ 3 ^ 2 for Xd <C x <C x coo i, where x coo i is the radius when the blast wave 
becomes strongly radiative (see §2.4). Asymptotes for N'(p';x) and U a di can easily be derived using this 
expression for P(x). Consider a blast wave that is initially relativistic so that Pq 3> 1. When x <C Xd, 
P(xi) ^ P , T(xi) = T(x) ^ r , and D = xjx,. One obtains 

N'(p'; x) = ^JV 2 x 3 H[p> - P ] , for x « x d , (21) 


and 

Uadi{x) = rn p irn x 3 P ai for x < x d . (22) 

The quantity H[a] = 1 for a > 1 and H{x) = otherwise. During the period preceeding deceleration, both 
the internal energy and total number of particles increases oc x 3 , and equation (21) shows that a quadratic 
tail of particles extending to lower energies is formed due to adiabatic losses. 

During the deceleration regime ^ « 1 <C T^ 3 Xd while the blast wave is still relativistic, P(x) = T(x) = 
Po(x/x d )~ 3 / 2 , and D = (x/x l ) 3 ^ 2 . One obtains 

N'(p'- x) - ^-5[p' P (-) 3/2 ] , for x d « x « T 2 Q /3 x d , (23) 

O Xd 
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and 

U adi {x) = {xx d fl\ for x d « x « T 2 /3 x d . (24) 

Equation (23) shows that during this phase, particles maintain a roughly monoenergetic distribution because 
adiabatic losses of particles balance the slowing down of the blast wave and the energy at which new particles 
are captured. At the end of the rclativistic deceleration phase at x = T a ' x d , the internal energy U ad i = E , 
that is, a large fraction of the fireball kinetic energy has been transformed into internal particle energy within 
the shell. 

The decelerating blast wave becomes nonrelativistic when x >• T^ 3 x d . In this regime, P{x) = 
Po(x/x d )~ 3 / 2 , T(x) = T(xi) = 1, and D = x/xi. The comoving distribution function and internal en- 
ergy of particles that are swept in after the blast wave becomes nonrelativistic are 

N'( P ';x) £ Svrno H \p; P (f )" 3/2 , P^fy], for x > T^x d , (25) 

and 

U adi (x) = nm p no P 2 x 3 [l - P Q 4/3 (^) 2 ], for x > T 2 /3 x d , (26) 

respectively. In equation (23), H[y; a, b] = 1 for a < y < b and H[y; a, b] = otherwise. Comparing equations 
(24) and (26), one sees that the internal energy approaches a constant value. 

The inset in Fig. 1 shows the numerical calculation of the internal energy U ad i- The above asymptotes can 
be seen to agree with these calculations. The difference between the adiabatic and momentum-conservation 
solutions can be appreciated by noting that in the absence of adiabatic losses, 

f* 

U mom {x) = A-Km p n / dxixf [T(xi) - 1] . (27) 
Jo 

Because particles experience no adiabatic losses in the momentum-conservation solution, the internal energy 
will always be greater than in the adiabatic case where such losses are included. Using equation (27), one sees 
that Uraom — (4/3)J7 a di when x <§C x d , and U mom = 2U ad i when x 3> x d noting that very little additional 
internal energy is swept into the blast wave when x > Iq S x d . During the episode of deceleration while 
the blast wave is relativistic, most of the internal energy is swept-up when X — 7 SO that p(xi) = P{x) 
and 7 = T(x) in equation (17). Thus the relativistic forms of the adiabatic and momentum-conservation 
solutions are similar. Most of the internal energy that remains within the blast wave as it decelerates to 
nonrelativistic speeds was introduced during the relativistic phase. The greater amount of internal energy 
for the momentum-conservation solution means that the blast wave must travel more slowly than in the 
adiabatic solution in order to conserve total momentum. 

The dependences in the various limits are obtained by noting that for the momentum-conservation 
solution, P{M + f* dx[dm(x)/dx]T(x)} ^ /3T[M + m(x)T(x)] ^ /3T[M + kx 3 T] = const, giving the 
asymptotes T oc x~ 3 / 2 when ro > T > 1 and (3 oc x~ 3 when r — 1 <C 1. For the radiative solution, 
P{M + f* dx[dm(x)/dx}} ^ PT[M + m(x)] = (3T[M + kx 3 ] = const, giving the asymptotes r oc x when 
r 3> T ^> 1 and (3 oc x~ 3 when r — 1 <C 1. The limits for the adiabatic solution would seem to be obtained 
through total energy conservation from the expression T{Mq+ Jq dxT(x)[dm(x) / dx}} = T[Mq + Tm(x)] = 
r[M + kx 3 (T — 1)] = const. This gives the correct relativistic asymptote T oc x~ 3 / 2 when r B » T » 1, 
but implies that (3 oc x~ 3 when r — 1 <C 1. As is apparent from the numerical results, the change in internal 
energy due to adiabatic losses becomes important in the nonrelativistic regime so that this estimate is not 
valid there. 
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2.4. Nonrelativistic Limit of Equation for Blast- Wave Evolution 

The evolution of the speed of a nonrelativistic blast wave is described by the equation 

d(Mp) 1 dU adi d/3 

-^r--p— + u Tf (28) 

which is obtained by letting P — > /3, T — > 1, and dx — > flcdt in equation (5). Because dU a< n/dt ~ (f3c)~ 1 U/x 
and fi <C 1, this term dominates the Udj3/dt term in equation (28). When consideration of the external 
medium pressure on blast- wave evolution is taken into account, we recover the equation of Cowsik and Wilson 
(1975). This was used to establish the important result that if SNe are the sources of cosmic rays below 
the knee of the CR spectrum, then CR acceleration must persist during the late phases of SNR evolution in 
order to overcome adiabatic energy losses within the expanding remnant. 

The dotted curve in Fig. 1 shows the dependence 

(Pa , for x < x d 

P{x) = ° = { (29) 

y/\ + {x/x A )* [ /3o( ^)-3/2 ; for Xd <<; x <<; Xcqo1 

suggested by the solution (Chiang and Dcrmcr 1999) to the equation for momentum conservation. This 
provides a representation that is accurate to within 10% of the numerical solution and the hydrodynamical 
asymptotes describing the behavior of an adiabatic blast wave. The Sedov length 

£ S ee T 2 /3 x d = (_?i?5_)V3 = 5.4 x 10 18 (^i) 1 ^ cm = 1.76(^i) 1 / 3 pc (30) 
47moTO p n n 

(Sari et al. 1996). 

At nonrelativistic energies, (Mo + 4irx 3 nom/3)v = MqVq from momentum conservation, where v — v(x) 
and Vo is the initial speed of the ejecta. Thus the SN ejecta begins to undergo significant deceleration when 
the swept-up mass is about equal to the ejecta mass, that is, when 47r£|nom = 3Mo- For To — 1 <C 1, this 
occurs when is ~ %d, recalling that E — > M in this limit. 

In the nonrelativistic limit with P <C 1, equation (29) reduces to 

{Po = v /c , for x < x d = £ s 

(31) 
/3o(f)- 3/2 = (|^) 1/2 ^ 3/2 , ior£ s «x«x cool 

where Ek e = 0qEq/2 is the kinetic energy of the injection event when (3q <§; 1, and po = mno is the rest-mass 
energy density of the CBM. From the nonrelativistic branch in equation (31), we see that 

p {x) ^0. 10 (^)^[x( P c)}- 3 / 2 , (32) 
no 

where E 5 i is the explosion kinetic energy in units of 10 51 ergs. For comparison, Lozinskaya (1992) gives 
the relation (3{x) = [4 x 2.02 x E ke /(25 x po)} 1/2 x-'^ 2 = 0.086(S 5 i/n ) 1/2 [a;(pc)]- 3/2 , obtained from the 
hydrodynamic self-similar Sedov solution. When x » x coo i, the Sedov phase ends. Chevalier (1974) and Fallc 
(1981) find that due to the onset of atomic and bremsstrahlung processes in the cooling plasma, x coo /(pc) 
= 19^i 29 n - 41 and x coo/ (pc) = 20££ 1 295 n ( 7°- 409 , respectively (see also Lozinskaya (1992)). 
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3. Particle Acceleration in GRB Blast Waves 

Particle acceleration can occur in GRB blast waves through a first-order Fermi mechanism involving 
internal or external shocks, and through second-order Fermi acceleration involving gyroresonant scattering of 
particles by magnetic turbulence in the magnetic field of the blast wave. Waxman (1995) and Waxman and 
Bahcall (2000) argue that shocks within a relativistic wind that persists over a timescale characteristic of the 
duration of a GRB could accelerate particles to UHECR energies. Such a blast wave must be highly radiative 
in order to give good efficiency for UHECR production; otherwise the bulk of the energy of the explosion 
is retained within the blast wave to be radiated during the afterglow phase when acceleration ceases and 
adiabatic losses reduce the energy of the nonthermal particles. UHECR acceleration in an internal shock 
model must also compete with strong photomeson losses that could limit acceleration to the highest energies. 

An external shock acceleration model for UHECR acceleration of particles accelerated from nonrela- 
tivistic energies has been shown to be infeasible (Gallant and Achterberg 1999). Following the first shock 
crossing, subsequent energy gains are modest because cosmic rays that diffuse in the surrounding medium 
are captured before pitch angle scattering changes the direction of the CR with respect to the shock normal 
by more than <~ 1/r. 

We therefore consider stochastic acceleration within the GRB blast wave. Wave turbulence in the 
magnetic field of the blast wave will be generated through the process of capturing electrons, protons, and 
charged dust (Pohl and Schlickeiser 2000; Schlickeiser and Dermer 2000). Rayleigh- Taylor instabilities in the 
expanding fluid and density irregularities in the surrounding medium will also introduce magnctohydrody- 
namic turbulence on sizes corresponding to the blast-wave width (Meszaros, Rees, and Papathanassiou 1994; 
Ryu and Vishniac 1987). This energy will cascade to MHD turbulence on smaller size scales. Although the 
wave generation and cascading process is inherently dynamic, we approximate the turbulence spectrum by 
a power-law in wavenumber space in order to provide a simplified treatment of the process. 

3.1. Maximum Energy of Particles Retained in GRB Blast Waves 

Allowed sites of UHECR acceleration must satisfy the Hillas (1984) condition, which essentially requires 
that the particle Larmor radius n, be less than the size scale D of the system. For GRB blast waves, 
tl = (m p /eB)p'(A/Z), where B is the mean magnetic field, Am p is the ionic mass and Ze its charge, and 
D = A'. When generalized for bulk relativistic motion of the acceleration region (e.g., Waxman (1995); 
Norman et al. (1995); Rachen and Mesazaros (1998)), the maximum measured particle energy that can be 
accelerated in a GRB blast wave is E max (ergs) = TeZBA' . 

The value of B is conventionally assigned in terms of a magnetic field parameter es that gives the 
magnetic field energy density in terms of the energy density of the downstream shocked fluid. Thus 

S(Gauss) = (327r^n e B m p ) 1/ V r ( r ~ 1) = 0.39(e B /m ) 1/2 iVry(l + r ) - 0.39(e B ^n ) 1/2 P , (33) 

where the approximation is valid to within a factor of \/2 from the extreme relativistic to the nonrelativistic 
regime. Recalling the definitions of A' in equation (11) and Xd from equation (19), we therefore see that 
p = Tp' < 0.39e(e By (mo) 1/2 ^/A2y r ( r - Tj/Am p , implying 

E max (eV) - 7.6 x 10 20 ^( T ^)(^2) 1/6 4 /2 (^54r 300 ) 1/3 (^)(-)v / l - • (34) 

1/lZ In Xd 
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The value of E max depends on x and P, but we can see that at the deceleration timescale 

*-3&- ,35) 

where T = To and x = xj, energetic GRBs can in principle accelerate protons and ions to ^> 10 20 eV. The 
weak dependence of E max on n shows that UHECR acceleration is still possible for n «l cm~ 3 , provided 
that e_e approaches its equipartition value (Waxman 1995). Equation (34) agrees with the expression derived 
by Vietri (1998a, 1995). 

The Hillas condition when applied to electrons with Z — 1 also yields equation (34). Once ions begin 
escaping from the acceleration region, however, electrostatic forces may develop to modify the lepton escape 
rate. Moreover, synchrotron radiation losses are much more important for leptons, as we see quantitatively 
in Section 3.2.6, and these losses will strongly impede lepton acceleration to energies set by this condition. 

The relationship between x and observer time t = t^T is obtained by noting that dx = PTcdt. Assuming 
an adiabatic blast wave, we can use equation (29) to obtain 

t , for t <C 1 

(4r)!/4 , fo r 1 « r « r* /3 (36) 

Xd { (5r/2r ) 2 / 5 , forr»r^ /3 , 

from which it follows that 

( 1 , for r « 1 

£_ o* J ( 4r )-3/8 ; for ! « T « r 8/3 (37) 

F ° { (5r/2r )- 3 / 5 , forr»r« /3 

for t < 1 

for 1 « r « rf 3 ( 38 ) 
for r > r* /3 . 

Expressions (36)-(38) apply to initially relativistic or nonrelativistic blast waves, but in the latter case the 
middle branches do not obtain insofar as T s=a 1. The late-time (r 3> T^ 3 ) behavior of the nonrelativis- 
tic expressions persists until the blast wave becomes highly radiative. In the nonrelativistic regime, the 
blast wave radius depends upon observer time according to the relation x = (75£t)C 2 /167mor7i p ) 1 / 5 i 2 / 5 = 
0.32(i?5i/no) ' 2 £ ' 4 (yr) pc, in good agreement with the Sedov- Taylor self-similar solution (Lozinskaya 1992). 

Equations (34-38) show that UHECR production in GRB blast waves is in accord with the Hillas limit 
throughout the prompt and afterglow phase of a GRB for suitably energetic GRBs when eg <~ 1. The <~ r~ 1//8 
dependence during the afterglow phase implies a reduction in E max by about one order of magnitude during 
the relativistic phase of a GRB blast wave, assuming that all other parameters remain constant. During the 
nonrelativistic phase, this dependence steepens to E max oc r" 1 / 5 . Insofar as the bulk of the total energy 
liberated by fireball transients and GRBs is emitted by the relatively infrequent, very energetic explosions 
with E 54 S> 0.1 (Dermer 2000), GRB blast waves satisfy the Hillas condition and are viable acceleration 
sites for UHECR production, provided that es approaches unity (see below). 



and 




3.2. Stochastic Particle Acceleration 



Our treatment of stochastic gyroresonant particle acceleration is very simplified. We make use of sys- 
tematic energy gain rates and diffusive escape timescales based upon expressions derived in the quasilinear 
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regime (Melrose 1974; Miller and Ramaty 1989; Bogdan et al. 1991), and consider only parallel-propagating 
modes that resonate with particles through the Doppler or anomalous Doppler resonance. The most impor- 
tant modes for protons and ions are shear Alfven waves, including ion-cyclotron waves, whereas electrons 
will gyroresonate with fast mode waves, including whistlers. We assume that the spectrum of parallel waves 
can be described by the function w(k) oc k~ q , where w(k)dk is the differential energy density in waves 
with wave numbers between k and k + dk and q is the spectral index of the turbulence. For a Kolmogorov 
spectrum of turbulence, (7 = 5/3. Assuming symmetry between the direction of wave propagation so that 
w(k) — w(—k), the ratio of energy density in waves to the magnetic-field energy density ub — B 2 /8tt is 
given by £ = 2 in dk w(k)/us- The smallest wave number k m i n is related to the largest size scale on 
which turbulence is generated. For a blast wave, we therefore assign k m i n w 1/A'. 

In the quasilinear regime, the systematic energy-gain rate for stochastic particle acceleration implied by 
the pitch-angle-averaged momentum diffusion coefficient with shear Alfven turbulence is given by 

(^)sto = |( i ^)/3iC(cfc mm )Kfc mm )"-V^ 1 (39) 



(Melrose 1974; Miller and Roberts 1995; Dermer et al. 1996), where c(3a = B / ^/Aim'm p is the Alfven speed 
and r° L = (m p /eB)(A/Z) is the nonrelativistic Larmor radius of the particle, so that n, = r®p'. Expression 
(39) also holds for gyroresonance with fast-mode waves away from the whistler branch. The magnetic-field 
prescription (33) implies that /3a — 2 3 / 2 e^ , 2 P. In order to accelerate UHECRs in GRB blast waves [eq.(34)], 
we require that es w 1, implying that /3a ~> 1- We interpret this unphysical result to indicate that the 
Alfven velocity approaches c whenever P > 1 and es ~ 1- We also assume that ( > 0.1, so that the magnetic 
field approaches the fully turbulent regime. When these conditions hold, stochastic particle acceleration is 
very rapid because the scattering centers are traveling at speeds approaching c, and the scattering rate is 
large. This is contrary to the usual situation found in the interstellar medium, where Pa ~ 10~ 4 and ( may 
be < 1. 

These conditions violate, however, the quasilinear assumptions used to derive the stochastic transport 
coefficients leading to equation (39). We nevertheless assume that this expression still provides the correct 
functional dependence on p' and q, so that 

p> sto = K p (c/A>)(r° L /Ay-yi-i , (40) 

where K p = Tr(q—l)j A /3 A (/3q generalizes equation (39) for gyroresonant interactions with relativistic plasma 
waves. The kinematic factor j A = 1/(1 — /3 A ) in K p allows for the possibility that the fractional change in 
momentum during the pitch angle isotropization time scale can exceed unity for relativistic Alfven waves. 
Consequently 

(dp , _ Psto _ / r LP \q-l ( A i\ 

1 dx )st ° - Pc~ Pr° L [ A' ' ■ [ ' 

Equation (40) implies an acceleration time scale 

tacc ~ l p> 1 ~ K p c [ A>> ~ K p c' (42j 

where the inequality holds because of the condition tl < A'. When K p — (27r)~ 1 , we recover the form 
that Rachen and Mesazaros (1998) use to gives the minimum acceleration time scale resulting from Fermi 
processes. Our subsequent comparisons agree with their conclusions when K p < (27r) _1 , in which case 
T 3> 10 2 , n > 1 cm~ 3 and es ~ 1 arc required to accelerate UHECRs. Stochastic gyroresonant acceleration 
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with relativistic Alfvcn waves can in principle permit K p to be 3> 1, in which case such extreme values of 
T, no, and es are not required. A detailed study of the resonant interactions of particles with relativistic 
plasma waves in a fully turbulent plasma is required to determine the range of possible values of K p , and is 
beyond the scope of the present study. Here we assume that K p ~ 1 in GRB blast waves. 

Particles will diffuse via pitch-angle scattering and subsequently escape from the blast-wave shell. The 
diffusive escape timescale along the direction of the large scale magnetic field is given in the quasilinear 
regime by 

t' esc « max[4,„, |(g - 1)(2 - g)(4 - q)t' dyn Q{r° L k min y- 2 p'«- 2 ] (43) 

(Barbosa 1979; Steinacker and Miller 1992; Schlickeiser 1989; Dermer et al. 1996), where t' dyn w A'/c 
represents the transit timescale for straight-line escape from the blast wave. Because much of the application 
in this paper falls outside the quasilinear regime, we again assume that the quasilinear proportionality holds 
and write 

t' esc *t' dyn m^[l,K t { r -£y- 2 ]. (44) 

The values that K t can assume are poorly known because K t depends on the magnetic field direction in 
the blast wave. If the magnetic field is primarily radial, then Kt could be -C 1, depending on the density 
of scatterers and the parallel diffusion coefficient. If there is a significant transverse magnetic field, then 
diffusive escape from the blast wave will be inhibited so that K t 1. 



3.2.1. Comparison between Available Time and Stochastic Gain Time Scales 

The highest energy a particle can reach is limited by the time available for particle acceleration. This 
maximum energy can be determined by solving the equation dp' /dx = (dp' /dx) st0l using equation (41) and 
noting that r^Y is independent of x for r > 1. Considering only proton (A/Z = 1) acceleration in the 
relativistic regime T ^> 1, we find that 

fc. a vT , + *,(MX10-).^(W^ _ , (45) 

J A 

where p\ and Xi refer to the injection momenta and location. For a Kolmogorov spectrum with q = 5/3, 
Pmax.av — 1-8 x 1 ~ 5 K p ^/cbJ^M) x . For standard parameters with x d = 2.6 x 10 16 cm and r = 300, 
Pmax,av — Fp' maxav — 1-4 x ^Q 15 K p \fe~B at x = x d . Thus the available time constraint does not prevent 
acceleration to ultra-high energies unless if p <l. 



3.2.2. Comparison between Adiabatic Loss and Stochastic Gain Rates 

Particles will be stochastically accelerated until adiabatic losses prevent further acceleration. This energy 
is given by the condition (dp'/dx) sto < \(dp' / dx) ad i \ = p'/x, where the last expression follows from equation 
(10), noting that \(dp' / dx) a di\ = (1 +g/3)/ x for T(x) cx x~ 9 . Even for a radiative blast wave with g = 3, the 
expression —p'/x for the adiabatic energy-loss rate is valid to within a factor of 2. This relation also roughly 
compares the timescale for a particle to be accelerated to some value of p' against the available time. One 
finds that the maximum momentum that can be reached before acceleration is halted by adiabatic losses is 

FA' 

p maxM i = Tp' max .a dl = (l2K p fl^ (-g-) . (46) 

'l 
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Recalling the Hillas condition p < TA'/rl, we therefore see that the maximum energy that a particle can 
attain, subject to adiabatic losses, is 

E m ax,adi (12-Kp) ^ ^ E max , (47) 

where we have set /a = 1/12. Thus we see that particles can be accelerated to energies given by equation 
(34) if K p > 1/12. It is also interesting to note that if q — 2, adiabatic losses will never dominate stochastic 
particle acceleration when K p > 1/12, and will always dominate acceleration if K p < 1/12. 

Equation (47) is in agreement with the results of Rachen and Mesazaros (1998) when allowance is made 
for the different approaches. Here we explicitly account for the spectrum of turbulence and determine A' 
from hydrodynamic considerations. Rachen and Mesazaros (1998) determine the size scale of the system 
from temporal variability, and define the adiabatic energy loss time scale in terms of \B/B\ which, if B is 
proportion to some power of x, yields equation (12) within a factor of order unity. Although our results are 
similar, our conclusions differ because, as noted above, we assume that K p can be > (2tt)^ 1 . 



3.2.3. Comparison between Diffusive Escape and Comoving Timescales 

When the timescale for a particle with a given momentum to escape diffusively from the blast wave 
is shorter than the comoving timescale, then acceleration to larger momenta is no longer possible. This 
condition is K t (A' /c)(r a L p' / A')"- 2 <t' = x/(cP), implying that Pmax . esc = (rA'/r£)(^ t / A ) 1 /(2-9). Conse- 
quently 

E max . esc > {K t l\lfl^E max . (48) 

The inequality applies when the transit timescale exceeds the diffusive escape timescale (eq. [44]). Equation 
(48) shows that if K t > Z^ 1 = 12, as could occur for large transverse magnetic fields in the blast wave, then 
particles will not diffuse from the blast wave during the available time before being accelerated to energies 
limited by the Hillas condition. 



3.2. 4- Comparison between Stochastic Energy Gain and Diffusive Escape Timescales 

Acceleration to a given energy will be limited by the diffusive escape from the blast wave. The inverse 
of the timescale for stochastic acceleration is t~^ Q = K p (c/ A')(r®p' / A') q ~ 2 . Comparing with the diffusive 
escape timescale (44) gives p < (IW /rl)(K p K t ) 1 ^ 4 - 2q \ or 

E i = IK R r ,1 1 /( 4 - 2 9) E (491 

^mai.occ/esc — \ x *-p*^t) '-'max ■ \^ I 

Thus we see that particles can be accelerated through stochastic gyroresonance to ultra-high energies if 
K p K t > 1, limited only by the condition that the Larmor radius is smaller than the blast wave width, if 
K p > /a and K t > 1//a- When q = 2, particles can either be accelerated to arbitrarily high energies if 
K p K t ^ 1, or will not be accelerated if K p K t < 1. 
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3.2.5. Comparison between Stochastic Energy Gain and Synchrotron Radiation Losses 

The synchrotron loss rate for particles with randomly oriented pitch angles in a magnetic field with 
mean strength B is given by -j' syn = 4co-tZ 4 (B 2 /8irm e c 2 )p' 2 /[3A 3 (m p /m e ) 3 }. Hence 



' S J/™ Q A3 f™ /™ ^2 ^ 1 syn ^ > 



16 /uZ 4 ccTTes/unor(r — 1) /2 
3 A 3 (m p /m e ) 2 

Equating this rate to the stochastic energy gain rate (40) gives the maximum particle momentum achievable 
due to the competition with synchrotron losses. The result is 

„ ^ v \ K p ( r l\«-2 3 ^ 3 (m p /m e ) 2 11/(3-9) 

Pma X ,syn-Ll A , ( A ,) W ^ a T e B fin T(T-l) 1 ' V L ) 

Figs. 2a-d show calculations of E max for protons implied by the Hillas condition, equation (34), compared 
with the maximum energy E maXySyn {eV) = 9.4 x 10 8 ' Ap. maXySyn determined by comparing the stochastic 
particle acceleration energy-gain rate with the synchrotron energy-loss rate from equation (51). Hence 
A = 1, Z = 1, and we let /i = 1. The standard parameter set employs a Kolmogorov turbulence spectrum 
with q = 5/3, cb = I^oo = ni = E^ = K p = 1, and /a = 1/12. The lower set of curves in each figure 
gives the Hillas condition result, and these curves are independent of the level of turbulence and particle 
acceleration rate set by K p . The limits set by available time, adiabatic losses, diffusive escape, and the 
particle acceleration rate give maximum energies proportional to E max , as seen in Sections 3.2.1 - 3.2.4. The 
upper set of curves gives the maximum energy E maXiSyn determined by the competition between particle 
acceleration and synchrotron losses. When these latter curves are less than E maxi then the maximum possible 
particle energy is set by E max , syn . 

For the standard parameter set, Xd = 2.6 x 10 16 cm and = ^o^ x d — 10 18 cm. The dissipation of 
directed kinetic outflow into internal particle energy is greatest for an adiabatic blast wave over this range 
of distances. Thus if GRBs are the sources of UHECRs, then the values of E max reached between Xd and is 
determine whether GRBs are viable sites of UHECRs. As can be seen, there are wide ranges of parameter 
values that permit particle acceleration to proton energies > 10 20 cV. Because more energetic ions can 
have comparable Larmor radii than less energetic protons due to their larger charges, ions can reach even 
larger energies than protons subject to the Hillas condition. The stronger synchrotron losses for ions are not 
sufficient to restrict the maximum particle energy of ions to values less than that for protons, because the 
gyroresonant waves are more effective in accelerating the ions, which have smaller Larmor radii for a given 
energy. One finds that E max . syn oc A^/^ (A/Z)^/^ . 

Competition of Fermi acceleration with synchrotron losses has been considered previously (Vietri 1995; 
Waxman 1995; Rachen and Mesazaros 1998). Other radiative losses, in particular, photomeson production 
of high energy particles interacting with synchrotron photons radiated by relativistic electrons in a GRB 
blast wave, can also prevent particle acceleration to ultra-high energies. Waxman (1995) and Rachen and 
Mesazaros (1998) consider neutrino production losses for general Fermi acceleration processes, Waxman and 
Bahcall (1997) consider neutrino production losses in a colliding shell model, and Vietri (1998b) considers 
neutrino production in a scenario where UHECRs are accelerated by first-order Fermi acceleration in rela- 
tivistic GRB shocks. However, Vietri (1998a) errs in his treatment of relativistic shock Fermi acceleration at 
an external shock (Gallant and Achterberg 1999), and so overestimates the highest energy neutrinos formed 
(Vietri 1998a, b). Photopion processes have been treated by Dermer (2000) in detail for an external shock 
model where UHECRs are accelerated through stochastic processes in the blast-wave shell, and shows that 
competition with photomeson losses does not prevent acceleration to ultra-high energies. 
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Figs. 2a-2c show that E max generally increases for increasing values of Tq, and 712, with all other 
parameters remaining constant. For the largest values of _E 54 and 712 that are considered, synchrotron losses 
rather than the Hillas condition determine the maximum particle energies. Fig. 2d illustrates the effects of 
changing K p , q, and es on E max . Changing K p and q has no effect on E max but does change E max ^ syn . 
If K p < 0.1, UHECR particle acceleration becomes infcasiblc in GRB blast waves. The large value of K p 
means that the plasma is fully turbulent and has relativistic Alfven speeds. Changes in E maXtSyn are only 
weakly dependent on the turbulence index q. By comparing the q = 5/3 and q = 2 results in equation (51), 
one finds that the parameter dependence typically varies oc 1/4 power. Because the q = 2 result is much less 
complicated, we examine the q = 2 results for Pmax,syn from equation (51). When r> 1, 

p m „ = 4-4 x 10 12 *E (-^-)- 1 A f 4 (^%-) 1/3 • (52) 

It is interesting to note that the standard parameters r 30 o = -E54 = n 2 = 1 used here are the same as 
those found to give good spectral fits to gamma-ray bursts (Dermer et al. 2000), and these parameters also 
give protons with energies 3> 10 20 eV. An important difference between parameters is that es — 10~ 4 in the 
spectral modeling, whereas here we let es — 1- The reason that es had to be so small in the spectral modeling 
was to avoid forming photon spectra resulting from cooling electron distributions, which are rarely observed 
in GRB spectra. But electrons will also be accelerated via gyroresonant stochastic particle processes. When 
particle acceleration is included, es can be much larger because acceleration competes with radiative losses 
and the formation of cooling electron distributions is avoided. 

In stochastic acceleration of highly relativistic electrons, helicity effects and gyroresonance with different 
plasma waves are not so important as for lower energy electrons where interactions with whistler waves are 
important (e.g., Levinson 1992), so we can apply equation (52) to determine the peak of the synchrotron 
spectrum radiated by electrons. Synchrotron losses rather than the Hillas condition determines the maximum 
energy of electrons accelerated through stochastic gyroresonance acceleration with plasma turbulence. The 
observed peak frequency of the vF v spectrum from accelerated electrons in a GRB source at redshift z is 
given by 

v max (Hz) = 2 A_^L BT{ P ^f (53) 

1 + z r 

for q = 2, where \f max is (m p /m e ) 2 smaller than p ma x,syn given by equation (52), and A = Z = 1. Thus 

tk t i^ 2 r 4 / 3 

E pk (keV) = fcw = _ (-) (_) 4/>/6nW ^ 

Equation (54) shows that E p k is in the range of E p k values observed with BATSE (Mallozzi et al. 1995) 
between w 100 keV and 1 MeV when r <~ 300-3000. This effect may resolve a fine-tuning issue in the 
external shock scenario (Chiang and Dermer 1999) whereby the value of es has to be large enough to 
provide good radiative efficiency, but not so large that cooling electron distributions are formed. 

Multiwavelength afterglow modeling within the standard blast wave model (Wijers and Galama 1999; 
Panaitescu and Kumar 2001) imply that n and e# are much smaller than the standard values used here, 
as already noted in Section 2.2. If these implied values are correct, then acceleration of UHECRs in GRB 
blast waves is very difficult during the afterglow phase, as is apparent from equation (34). For example, 
consider the parameters derived by Wijers and Galama (1999) for GRB 990508, namely no = 0.03 cm~ 3 , 
E 54 — 0.03, and cb = 0.12. Acceleration of > 10 20 eV protons requires unreasonably large initial Lorentz 
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factors, namely r ~ 3 x 10 4 . Similarly conclusions are obtained using parameters for GRB 980703 and 
GRB 990123 derived by Panaitescu and Kumar (2001). 

This may imply that the bulk of UHECR acceleration is done by GRBs during the prompt phase where 
the derived parameters may not hold. The inclusion of stochastic particle acceleration within the blast wave 
itself represents, moreover, a fundamental departure from the standard blast wave model. Acceleration and 
radiation losses take place concurrently in the model under consideration, as compared with the standard 
assumption where particles are instantaneously accelerated and subsequently cool. Moreover, escape of 
UHECRs from the blast wave may drive the system into different radiative regimes that will complicate 
parameter estimation. An important test to check the derived value of the blast-wave magnetic field is to 
observe the appearance of the synchrotron-self Compton component in X-ray afterglow spectra, because its 
intensity is very sensitive to eg (Dcrmer et al. 2000). 



4. Solutions to Particle Continuity Equation 

We examine particle acceleration, escape, and losses within a GRB blast wave by employing a particle 
continuity equation 

^ + |^~ + SI^'>' <55) 

The term p represents the total energy-change rate of particles due to acceleration, adiabatic and radiative 
processes, and t esc is the particle escape timescale from the blast wave. Primes have been dropped to simplify 
the notation. The analytic solution to this equation when p = p(p) < and t esc — t esc (p) are independent 
of t is 

»( f ;„^[ W , 1 >p[-f^l (50, 

where t* — t — JJ dp" /\p(p")\. The solution for p > is obtained by reversing the limits on the integrals in 
the exponential and in the definition of t* , and by integrating from to p in the integration over p* . 

When the terms p and t esc depend on both p and t, the solution to equation (55) can be obtained 
semi-analytically by inspection, giving 

»<* " " jf * * (p ' *>' «*" jf u.Mp"t,r),t'\ } ■ <57) 

where pi = Pi(p, t, t") is obtained by inverting p = p(pi, t, t"), which solves dp/dt = p. Equation (57) can be 
shown to reduce to equation (56) when t esc and p are time-independent. 

For particle acceleration in a blast wave, it is more convenient to consider the evolution of the comoving 
distribution function with location x. Using the relationship dx = Pcdt' between radius and comoving time, 
equation (55) becomes 

dN(p;x) d p m dp N(p;x) 

~^x— + dp- {[ p- c + (^-d^te*)) + Pct esc ( P ,x) = Q{p ' x) • (58) 

The term p m represents the sum of the particle energy gain and loss rates, excluding adiabatic losses which 
are given through the term (dp/dx) a di from equation (12). Following equation (57), the solution to equation 
(58) is 

f x dpi f x dx* 

N(p;x) = 47rno / dxiCix^Xi 6\pi(p,x,Xi)- P(xi)]\-r±\ exp{- — — r— T — . (59) 

Jo d P Jxi Pct esc [p t (p,x,x*),x*\ 
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Here we substitute equation (15) for the source function Q, though modulated by the factor C(xi). Equation 
(59) can be solved analytically in some restricted cases, but must be solved numerically for the general case 
involving adiabatic losses and stochastic energy gains. When the energy lost through particle escape is a 
large fraction of the swept-up energy, the system departs from its adiabatic behavior. In this case, a more 
general numerical approach involving inversion of a tri-diagonal matrix must be used. We describe the 
numerical considerations for dealing with this system, though here we only report the results for blast waves 
that evolve in the adiabatic regime. 



4.1. Analytic Solution 

Equation (59) can only be solved analytically when the turbulence index q = 2. Depending on how 
the MHD turbulence develops and evolves, this case may be of primary interest. It is generally thought 
that wave energy is introduced on the largest length scales of the system and is transferred to smaller 
scales through nonlinear mode coupling (e.g., Zhou and Matthaeus 1990), to be dissipated through particle 
acceleration when the particle Larmor radii are in resonance with cascading turbulence. In the Kolmogorov 
phenomenology, the wave energy achieves a steady-state turbulence spectrum with index q = 5/3. 

Particles can only be accelerated when (dp/dx) s to > \{dp/dx) a di\ (Section 3.2.2), that is, when p > 
(A/r£)(/A/^p) 1/(9_2) if q > 2, and when p < (A/r£)(/ A / Kp) 1 ^ 2 -^ if q < 2. If q > 2, most of the MHD 
turbulence energy is contained in small wavevector (long wavelength) turbulence that accelerates high-energy 
particles with large Larmor radii, and vice versa for q < 2. Consequently, there will be a resonance gap that 
prevents particle acceleration when q > 2, because adiabatic losses dominate particle acceleration except 
for the very highest energy particles, of which there may be none if (A/rj- )(/a/ Kp) 1 /^ 2 ^ > P . If the 
turbulence spectrum evolves to the Kolmogorov index through a progressive hardening of the turbulence 
spectrum, then the low-energy swept-up particles will begin to sample the evolving turbulence spectrum 
when q = 2. If K p > /a and q = 2, stochastic acceleration will dominate adiabatic losses and accelerate 
particles to the highest energies. The waves will be damped by transferring their energy to the particles. In 
this way, a quasi-steady turbulence spectrum with q = 2 may be formed. 

When q = 2 and K p > /a, the stochastic energy gain rate dominates the adiabatic energy- loss rate for 
all particle momenta. Thus (dp/dx) s to — K p p/(PA) from equation (41). In the limit P ^> 1, (dp/dx) s to — 
K p p/(f&x) and, provided K t > 1, t esc = K t A/c from equation (44). The particle momenta therefore evolve 
according to the relation p = p i (x/x i ) K " . Equation (59) becomes 



N(p;x) = 47m f dx { C(x 2 )x 2 S[p(^) K ^^ - P(x t )} (£-)-(*p+V** 
Jo x Xi 



(60) 



There is no high-energy cutoff to the accelerated particle momentum in this limit, because the acceleration 
rate (dp/dx) sto oc 1/x and thus rapidly accelerates particles at the earliest times. 

When i; < id, P{xi) = Po and 

N(p; x) = inno C{ Xi ) (^) x 3 (%+3A>/*p+V(*p*0 , (61) 
K-pPa p 

where x t = x(^-) f&/K ". When x { > x d , P(xi) = P (xi/x d y 3 / 2 and 

Attuq Cjxi) xj ( x / Xl )-(^/^)/U 
[P ' ] lK pP (x l /x)K P /f-/f A x l } + l3Po(i l /xy3/ 2 /2x l }' [ } 
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where Xl = {P x 3 d /2 x K »/f* / p )V(*p//a+3/2) 

Fig. 3 shows the comoving particle distribution function that results from stochastic gyroresonant ac- 
celeration by plasma turbulence with q = 2 and K p — 0.5. We let K t = 2 in Fig. 3a and K t — 10 in Fig. 
3b. The other parameters are the same as used in Figs. 1 and 2, namely E54 = ni = = e# = 1, and 
we let the modulation factor C(xi) = 1. The particle distribution is multiplied by the factor m e c 2 pj in 
order to indicate the energy contained in the accelerated particles. The comoving particle distribution is 
well-described by a power law spectrum, and the index depends sensitively on the assumed values of K p 
and K t - For these parameters, particles can reach ultra- high energies before the blast wave decelerates to 
nonrelativistic speeds, which occurs at x = 45xd for the chosen parameters. Because the escape timescalc 
is independent of particle energy, power-law distributions of particles also escape from the blast wave. The 
highest energy escaping particles with large Larmor radii can freely escape to become UHECRs. In contrast, 
escaping particles with lower energies and smaller Larmor radii will be subject to the cosmic-ray adiabatic 
loss problem if their energy density exceeds the magnetic-field energy density in the local ISM. 



4.2. Numerical Solution 

To solve the continuity equation (59) numerically we discretized this equation using a conservative finite 
differencing scheme. To this end we introduce 



Equation (59) becomes 



ATi+l ATI + 1 _ ATi+1 



N) = N(xi, Pj ) ; F) = p{xi,pj)N(xi,pj) . (63) 
r i + 1 jyi+l 

1^111 til + Qi (64) 

Ax Ap t esc (xi,pj) 3 

We make the identifications N* +1 — Nf and N* +1 — N? . This leads to larger diagonal elements in the 
resulting system of equations and enhances the stability of the numerical algorithm. As a result, we are led 
at each x-step to a tridiagonal system of equations for AH, namely 

a 3 N^+b 3 N^ 1 + Cj N^\=S) 1 (65) 

where a 3 - = 0, bj ; = 1 + Ax[t~} c { Xil pj) — p/ Ap], Cj = pAx/Ap, and 5* = Q % - + N 1 -. In our computations, 
the grid in p consisted of 10 4 intervals, with a linear spacing in the range 1 < p < 10 5 , and a logarithmic 
spacing in the range 10 5 < p < 10 9 . The step size in x Ax = xj + i — X j was determined adaptively at each 
step to ensure that system (65) is well-behaved. The resulting Ax varied between 10 10 and 10 12 cm during 
the simulations. 

Figs. 4 and 5 show the numerical simulation results for a Kolmogorov turbulence spectrum q = 5/3 
when the blast wave reaches different locations between O.lxd and lOx^. The evolution of the blast wave is 
assumed to follow the adiabatic behavior described by equation (29). We use K p = 0.1, and K t — 10 in Fig. 
4, and K p = 1 and K t = 10 in Fig. 5. The other parameters arc the same as used in Fig. 4. Panel (a) in the 
figures show the time evolution of the comoving particle distribution, and panel (b) shows the cumulative 
energy distribution function of escaping particles measured by an outside observer, obtained by solving 

at ( \ r j- N esc \p esc /T(x),x] 

N esc (p esc ;x) = / dx — r , (66) 

Jo 1 \ x ) 

where N esc (p,x) — N(p; x)/ Pct esc (p, x) from equation (58). 
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In Fig. 4, we let the modulation factor C(xi) = 1. The development of the internal comoving particle 
distribution is shown in Fig. 4a at different blast-wave locations. The particles diffuse to higher energies with 
time, reaching a maximum value that is about an order-of-magnitude less than the maximum momentum 
permitted through the available-time constraint that was estimated from equation (45). This equation does 
not, however, take into account adiabatic losses which appreciably retard acceleration of the highest energy 
particles. Fig. 4b shows the escaping particle distribution that would be measured in an external system. 
Most of the energy is carried by the very highest energy particles with energies between 10 18 and 10 20 eV. 

The escaping UHECRs carry more than 10 56 ergs of energy, because wave damping and energy con- 
servation have not been taken into account in the calculation. A realistic calculation must consider these 
processes, as in treatments of impulsive Solar flares (e.g., Miller 1998; Miller and Roberts 1995). Moreover, 
we have restricted the blast wave evolution to the adiabatic regime. The calculation of Fig. 4 does not 
violate energy conservation or the assumption of an adiabatic blast wave if we assume that C{xi) < 10~ 3 , 
so that only a very small fraction of swept-up particles are accelerated to the highest energies. It is artificial 
to assume that C(xi) does not change with time location, but a self-consistent treatment of wave cascading 
and damping, required to obtain a correct evaluation of C(xi), is beyond the scope of this paper. 

In Fig. 5 we show a calculation of the evolving particle distribution with K p = 1. Here we assume 
that C(xi) = 10~ 4 so that the calculation is self-consistent. The effect of increasing the rate constant K p 
is to permit particle acceleration to much higher energies, until adiabatic losses and escape suppress further 
acceleration. Fig. 5b shows that a significant fraction of the energy is carried by particles with energies 
> 10 20 eV. In this calculation, the cumulative escaping particle spectrum develops a spectrum with number 
index <~ —3. Larger values of K t or other parameters of the calculation, such as r , can yield even higher 
energy escaping particles. Thus stochastic particle acceleration in relativistic blast waves can, in principle, 
produce cosmic rays with energies > 10 20 eV. 

4.3. Cosmic Ray and Ultra-High Energy Cosmic Ray Origin 

As shown in Figs. 4 and 5, protons accelerated to > 10 19 eV by a Kolmogorov MHD turbulence spectrum 
can diffusively escape from the blast wave during the prompt and afterglow phases of a GRB to become 
UHECRs. The Larmor radii of these particles are so large that they do not strongly couple with the 
interstellar gas and thereby avoid further adiabatic losses. This is because the energy density of UHECRs 
within a volume defined by their Larmor radius is smaller than the energy density of the mean magnetic 
field within that volume. UHECRs consequently behave as isolated particles rather than as a relativistic 
fluid that does work on its surroundings. 

This is not the case for lower energy particles which, if they escaped into the ISM, would generate 
streaming instabilities that would strongly couple these particles with the ISM. The adiabatic loss problem 
would be severe for these particles. As shown by the simulations, however, only the very highest energy 
particles escape during the relativistic phases of the blast wave. Once the blast wave enters the nonrelativistic 
Scdov phase, stochastic acceleration to ultra-high energies is much weaker because the Alfven speed becomes 
■C c. During this phase, shock Fermi acceleration becomes more efficient than stochastic acceleration. The 
internal particle population confined by the blast wave loses energy due to adiabatic expansion of the blast 
wave shell, but also begins to be energized through shock acceleration. The nonrelativistic phase of a GRB 
remnant evolution resembles a very energetic SNR, though differing from a conventional SNR by containing 
a population of high energy particles available from the earlier acceleration episode. Consequently the 
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timescale constraint (Lagage and Cesarsky 1983) on shock Fermi acceleration is relaxed. Cosmic rays will 
leave the remnant as the SNR dissipates in the ISM so that the cosmic-ray adiabatic loss problem is solved, 
as in the standard model, by shock acceleration taking place throughout the expansion of the SNR until the 
energy density of the relativistic particle fluid is small compared to the magnetic-field energy density of the 
ISM. 

Future work must treat both first- and second-order particle acceleration in a relativistic blast wave that 
decelerates to nonrelativistic speeds and determine whether this scenario can satisfactorily explain the origin 
of the hadronic cosmic rays above the knee of the cosmic ray spectrum, as well as making a contribution to 
hadronic CR origin at lower energies. 

5. Summary and Conclusions 

In this paper, we have treated adiabatic losses stochastic particle acceleration in blast waves formed 
by the ejecta from exploding stars. Because the Alfven speed approaches c whenever the magnetic field 
approaches its equipartition value in relativistic shocks, second-order Fermi acceleration can be much more 
efficient than first-order Fermi acceleration in GRB blast waves. In nonrelativistic SNR shock waves, by 
contrast, the first-order process dominates. Consequently, we have considered particle acceleration through 
stochastic gyroresonant acceleration in GRB blast waves during their relativistic phases. Even in the highly 
simplified approach developed here, it was necessary to treat adiabatic losses correctly in the expanding blast 
wave. This was dealt with in Section 2, where an equation for blast-wave evolution was derived (equation 
(5)) that contains a term that rechannels the energy lost by adiabatic expansion of the nonthermal particles 
into the directed outflow of the blast wave. An expression for the energy lost through adiabatic expansion 
by a general particle distribution function that contains both relativistic and nonrelativistic particles was 
derived in Section 2.2. The evolution of the internal particle distribution due to adiabatic losses was treated 
in Section 2.3, and an expression describing the evolution of an adiabatic blast wave which agreed with 
numerical calculations, was shown to reduce to well-known results in the nonrelativistic regime in Section 
2.4. 

In Section 3 we showed that ultra-high energy particle production in GRB blast waves satisfies the Hillas 
condition, which compares the particle Larmor radius with the size scale of the system. In the case of the 
GRB system, the size scale is the blast-wave width. Thus GRBs are allowed sites for UHECR production if 
the magnetic field approaches its equipartition value. We employed simplified forms for the energy gain rate 
and escape time scale due to stochastic gyroresonant interactions of particles with MHD turbulence based 
upon expressions derived in the quasilinear approximation. In Section 3.2, we derived the conditions under 
which competition between stochastic acceleration, available time, adiabatic losses, and diffusive escape 
permits acceleration to ultra-high energies. The conditions under which synchrotron losses limit particle 
acceleration were also examined. We find that wide ranges of parameter values permit particle acceleration 
to ultra-high energies through stochastic gyroresonant processes in GRB blast waves. However, if parameters 
derived in afterglow modeling are correct, then many GRBs cannot accelerate UHECRs during the afterglow 
phases. In Section 4, we numerically calculated comoving and escaping particle distributions by employing a 
continuity equation, and found that a population of UHECRs will diffusively escape from GRB blast waves 
during the prompt and afterglow phases. 

This study is preliminary and does not address the question of the origin of the turbulence, nor wave 
cascading and damping. We also assume that the blast wave is homogeneous, whereas it is in fact subject to 
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instabilities and nonlinear effects on the shock structure from the accelerated particle distribution. Within 
the constraints of the problem that are rather well defined, for example, the maximum blast-wave magnetic 
field, the characteristic blast-wave width, the adiabatic loss rate, and the time available for acceleration, 
we conclude that particle acceleration to ultra-high energies is possible in GRB blast waves. This process 
therefore offers a possible solution to UHECR acceleration and the origin of hadronic cosmic rays above the 
knee of the cosmic-ray spectrum. 

We thank A. Achterberg, M. Baring, L. O'C. Drury, A. Levinson, J. A. Miller, and H. Volk for dis- 
cussions, and the referee for constructive criticisms. The work of CD. is supported by the Office of Naval 
Research and the NASA Astrophysics Theory Program (DPR S-13756G). 
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Fig. 1. — Evolution of blast wave momentum P for different assumptions for the particle energy losses. 
The solid curve labeled "full treatment of adiabatic losses" is a numerical solution of equation (5) that 
includes adiabatic losses of swept-up particle energy that is transformed into the directed kinetic energy of 
the outflow. The curve labeled "momentum conservation solution" neglects adiabatic energy losses of the 
swept-up particles, and the curve labeled "radiative solution" assumes that the internal energy is promptly 
radiated. Dashed lines are hydrodynamic approximations at relativistic and nonrelativistic momenta, and 
the dotted curve is equation (29). The inset shows the evolution of the internal energy U mom and U a( n for 
the momentum-conservation and adiabatic solutions, respectively. 
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Fig. 2. — Maximum proton energies E max .ji as a function of blast wave location x implied by the Hillas 
condition, equation (34), obtained by comparing the Larmor radius and blast wave width, and E maXtSyn 
from equation (51), obtained by comparing the stochastic acceleration rate with the synchrotron cooling 
rate. Standard parameters are shown in the figure legend. Changes in E maXt n and E. maXySyn for different 
values of r , £54, and no are shown in (a), (b), and (c), respectively. Changes in E max , syn for different 
turbulence spectral indices q and acceleration rates K p , and &b are shown in (d), as well as changes in 
E m ax,n and E maXtSyn due to changing e B values. 
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Fig. 2.— Fig.2c (left), Fig. 2d (right) 
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Fig. 3. — Evolution of comoving particle distribution at different locations x of the blast wave in units of 
the deceleration radius Xd for the analytic result, equations (61) and (61), for a turbulence spectrum with 
q = 2. Parameters are -E54 = 1, no = 100 cm~ 3 , Tq = 300, and /a = 1/12. In panel (a), K p = 0.5 and 
K t = 2. In panel (b), K p = 0.5 and K t = 10. 
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Fig. 4. — Numerical calculation of the comoving particle distribution in panel (a), and the escaping particle 
distribution in panel (b) at different blast-wave locations. Parameters are the same as in Fig. 3, except that 
q = 5/3, K p = 0.1 and K t = 10. 
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Fig. 5. — Same as Fig. 4, except that K p = 1 and K t = 10. 



